


clear all
set more off
pause off
set matsize 11000



capture confirm file "Main_Results/"
if _rc mkdir "Main_Results/"

capture confirm file "Main_Results/Fig3/"
if _rc mkdir "Main_Results/Fig3/"

capture confirm file "Main_Results/Fig3/gphs/"
if _rc mkdir "Main_Results/Fig3/gphs/"

capture confirm file "Main_Results/Tables/"
if _rc mkdir "Main_Results/Tables/"


local region  Euro 

local response R_real




***************************************
***************************************
* LP from when the pandemic ended
***************************************
***************************************
use "Data/Big_Pandemics.dta", clear
gen pandemic=1
drop year
ren endyear year
sort year
merge 1:1 year using "Data/RealRateSince1310.dta"
sort year
replace pandemic=0 if pandemic==.
gen pandemic100k=cond(dead>=100000,pandemic,0)
cap drop Horizon
gen Horizon = _n-1 if _n<=41
merge 1:1 year using "Data/sdeaths.dta" , nogen


gen war_deaths = sdeaths
replace sdeaths =0 if sdeaths==.
replace sdeaths= cond(sdeaths>20000,1,0)
sort year

tsset year

cap drop b0* se0* 
cap drop h0* l0* h1* l1*
tsset year
foreach x of local region {
	foreach y of local response {
	
	gen b0_`x'`y'=0 if Horizon==0
	gen se0_`x'`y'=0 if Horizon==0	
	gen h0_`x'`y'=0 if Horizon==0
	gen l0_`x'`y'=0 if Horizon==0
	gen h1_`x'`y'=0 if Horizon==0
	gen l1_`x'`y'=0 if Horizon==0
	gen D_b0_`x'`y'=0 if Horizon==0
	gen D_se0_`x'`y'=0 if Horizon==0	
	gen D_h0_`x'`y'=0 if Horizon==0
	gen D_l0_`x'`y'=0 if Horizon==0
	gen D_h1_`x'`y'=0 if Horizon==0
	gen D_l1_`x'`y'=0 if Horizon==0
	gen pval_test_`x'`y' = 0 if Horizon==0
	
	cap drop raw
gen raw = .
replace raw = `x'`y'
/********** State space specification starts here *************/

constraint 1	[rstar]L.rstar		= 1		/* Random Walk */	
constraint 2	[raw]rstar			= 1		/* Definitional */


cap drop zero								/* generate zero line for figures */
gen zero = 0

tsset year
cap drop r_star

/************	 State space equations	 ******************/

sspace		(rstar L.rstar, state noconstant)		/// /* r-star */
			(raw rstar, noconstant ),				/// /* raw */
			covstate(identity) constraints(1/2) difficult /*method(kdiffuse)*/ iterate(25)
			
			qui{
			/* Generate smoothed estimates of r-star */
			predict r_star if e(sample), states smethod(smooth) equation(rstar)
			
	

			cap drop `x'`y'`trend
			gen	`x'`y'trend = r_star
			}

			
			


								}
								
						}


foreach y of local response {
	foreach x of local region {
	eststo clear
		forv h=0/40 {
		
						cap drop  d`y'`h'
						gen d`y'`h' = f`h'.`x'`y'trend- l.`x'`y'trend
						reg d`y'`h' pandemic100k sdeaths l(1/10).pandemic100k  l(1/10).sdeaths l(2/11).`x'`y'trend  if year>=1495
						
						eststo ols`h'
												
					
						replace b0_`x'`y'  = _b[pandemic100k]                      if Horizon==`h'
						replace se0_`x'`y' = _se[pandemic100k]                     if Horizon==`h'
						replace h0_`x'`y'  = _b[pandemic100k]+2*_se[pandemic100k]  if Horizon==`h'
						replace l0_`x'`y'  = _b[pandemic100k]-2*_se[pandemic100k]  if Horizon==`h'
						replace h1_`x'`y'  = _b[pandemic100k]+1*_se[pandemic100k]  if Horizon==`h'
						replace l1_`x'`y'  = _b[pandemic100k]-1*_se[pandemic100k]  if Horizon==`h'

						replace D_b0_`x'`y'  = _b[sdeaths]                      if Horizon==`h'
						replace D_se0_`x'`y' = _se[sdeaths]                     if Horizon==`h'
						replace D_h0_`x'`y'  = _b[sdeaths]+2*_se[sdeaths]  if Horizon==`h'
						replace D_l0_`x'`y'  = _b[sdeaths]-2*_se[sdeaths]  if Horizon==`h'
						replace D_h1_`x'`y'  = _b[sdeaths]+1*_se[sdeaths]  if Horizon==`h'
						replace D_l1_`x'`y'  = _b[sdeaths]-1*_se[sdeaths]  if Horizon==`h'
						
						
						

					
					
						test pandemic100k = sdeaths
						replace pval_test_`x'`y' = r(p) if Horizon==`h'
						estadd scalar pval_wars_pandemics = r(p): ols`h'
}
						
		esttab, keep(pandemic100k sdeaths)

		cap drop Zero
		gen Zero=0
		twoway ///
		(  rarea h0_`x'`y' l0_`x'`y' Horizon if Horizon~=. , fc(blue)  fi(10) lcolor(white)  lpattern(none)  )  ///
		(  rarea h1_`x'`y' l1_`x'`y' Horizon if Horizon~=. , fc(blue)  fi(20) lcolor(white) lw(none)  lpattern(none)  )  ///
		(  rarea D_h0_`x'`y' D_l0_`x'`y' Horizon if Horizon~=. , fc(red%50)  fi(10) lcolor(white)  lpattern(none)  )  ///
		(  rarea D_h1_`x'`y' D_l1_`x'`y' Horizon if Horizon~=. , fc(red%50)  fi(20) lcolor(white) lw(none)  lpattern(none)  )  ///
		(  line  Zero    Horizon if Horizon~=. , lc(gray)  lw(med) lp(dash)                     )  ///		
		(  line  D_b0_`x'`y'     Horizon if Horizon~=. , lc(red)  lw(thick) lp(dash_dot)  )  ///
		(  line  b0_`x'`y'     Horizon if Horizon~=. , lc(blue)  lw(thick) lp(solid)  ///
				xlab(0(10)40, labsize(large)) ylab(-2(1)2, angle(0) labsize(large)  nogrid) xtitle("Years since pandemic event ended", size(large)) ///
				ytitle("Percent", size(large))) ///
		, title()    ///
		xsize(5) ysize(4) scale(1.1) ysc(r(-2.5 2.1))  scheme(s1color)  ///aspect(0.5) scale(0.75)
		legend(row(2) order(6 7 ) region(color(none)) size(vlarge)  label(6 "wars") label(7 "pandemics") ring(0) position(11) ) ///
		name(`x'`y'_wars,replace)  
	
	
	
graph save   Main_Results/Fig3/gphs/Fig3b.gph, replace
graph export Main_Results/Fig3/Fig3b.pdf, replace

}
}






/* // If you want to export the coefficients into excel
keep Horizon b0* h0* h1* l0* l1* D_b0* D_h0* D_h1* D_l0* D_l1* pval_*
drop if Horizon==.
export excel  Horizon b0* h0* h1* l0* l1* D_b0* D_h0* D_h1* D_l0* D_l1* pval_* using "Main_Results/coefs_for_r-star", firstrow(varlabels) replace
*/

/* Output a table with coefficients */
do Programs/make_table

esttab pandemics wars pval_test_wp  ///
	using Main_Results/Tables/Table1a.tex, page r ///
	se(2) b(2) star(* 0.1 ** 0.05 *** 0.01) ///
	coef(h0 "$ h = 0 $" h10 "$ h = 10 $" h20 "$ h = 20 $" h30 "$ h = 30 $" h40 "$ h = 40 $" ) ///
	mtitles("Pandemics" "Wars" "p-value") ///
	stats(Obs, fmt(0) l("Observations"))  








